library(sfsmisc)
library(igraph)
alpha<-.5;min<-1;max<-10000;n<-100000
setwd("C:\\Users\\Ofer\\Desktop\\My Dropbox\\PhD\\LATEX\\4\\figures")

f.dist<-function(x, min=1, max=10, alpha=1.4){alpha*x^(-alpha-1)/(min^-alpha - max^-alpha)}
f.cum<-function(x, min=1, max=10, alpha=1.4){(x^-alpha-max^-alpha)/(min^-alpha - max^-alpha)}
f.cum.app<-function(x, min=1, max=10, alpha=1.4){x^-alpha*(min^alpha*max^alpha)/(max^alpha - min^alpha)}
generatePower<- function(nmax, alpha, min=1, max=Inf) {
	u<-runif(nmax) # alpha<-1; min=1; max=10; u<-.001
	norm<-(min^-alpha - max^-alpha)
	(min^-alpha - u*norm)^-(1/alpha)
}
d<-generatePower(n, alpha, min=min, max=max)
# hist<-hist(d, plot = TRUE, breaks=n/10, freq=F, xlim=c(1,20), add=T, border="black")
hist<-hist(d, plot = FALSE, breaks=n, freq=F)
normalized<-rev(cumsum(rev(hist$density)))/sum(hist$density)

########################################################
# exact function of the distribution
# png("ImgPowerDist.png",  bg = "transparent")
plot(seq(min, max, length.out=n), f.dist(seq(min, max, length.out=n), min=min, max=max, alpha=alpha), 
	# from=min, to=max, 
	col="black", lwd=2,  type='l',
    	xlim=c(0,20),
	ylab='',  lty = 2, 
	xlab='', pch=20, cex=.1,  xaxt = "n",  yaxt = "n")

mtext("frequency", side=2, line=2.5, cex=1.5)
mtext("x-value", side=1, line=2.5, cex=1.5)
eaxis(1, cex.axis=.8)
eaxis(2, cex.axis=.8)

# integrate(f.dist, lower=min, upper=max, min=min, max=max, alpha=alpha)
# plot distribution
points( hist$mids, hist$density, 	col="blue", pch=3, cex=1, lwd=1.5)
#legend
leg.txt <- c("Exact function (DF)",
		 "Simulation     (DF)")
legend(x=5, y=.5, leg.txt,  lty = c(2, -1), pch = c(-1, 3),
      col = c("black", "blue"), 
	cex = 1, lwd=c(2,1))

dev.off()

################################
plf<-power.law.fit(d, xmin = 1)
summary(plf)
?mle

#############################################################
# plot log cumulative distribution
# png("ImgPowerDistLogLog.png",  bg = "transparent")

curve(f.cum(x, min=min, max=max-1e-9, alpha=alpha), 
	from=min, to=max, col="red",
    	# add=T,
	ylim=c(1e-5, 10), lwd=2,
	ylab='', xlab='', 
	pch=20, cex=.01,  xaxt = "n",  yaxt = "n" , log="xy")
mtext("frequency", side=2, line=2.5, cex=1.5)
mtext("x-value", side=1, line=2.5, cex=1.5)
eaxis(1, cex.axis=.8)
eaxis(2, cex.axis=.8)


# plot log cumulative distribution
t0<-seq(log(1), log(length(hist$mids)), by=log(length(hist$mids))/100)
t<-as.numeric(as.vector(as.data.frame(table(floor(exp(t0))))[,1]))
#plot( hist$mids[t], normalized[t], log="xy", col="blue", pch=4, cex=1, ylab='', xlab='x value', xaxt = "n",  yaxt = "n")
#lm<-lm(log(normalized[1:60])~log(hist$mids[1:60]))
#summary(lm)


points( hist$mids[t], normalized[t], 
	col="blue", pch=4, cex=1, ylab='', xlab='x value', xaxt = "n",  yaxt = "n")
points(exp(hist$mids[t]-hist$mids[1]), exp(-0.5*(hist$mids[t]-hist$mids[1])), type="l")

# plot log distribution
points(hist$mids, hist$counts/n,  
	# xaxt = "n",  yaxt = "n", log="xy", 
	col="black", pch=3, cex=1)

curve(f.dist(x, min=min, max=max-1e-10, alpha=alpha), col="black", lty=2,
	pch=20, cex=.01, add=T, lwd=2)
#legend
leg.txt <- c("Exact function (DF)",
		 "Simulation     (DF)", 
		 "Exact function (reverse CDF)", 
             "Simulation     (reverse CDF)", 
		 "linear approx (reverse CDF)" )
legend(x=50, y=10, leg.txt, lty = c(2, -1, 1, -1, 1), pch = c(-1, 3, -1, 4, -1),
      col = c("black", "black", "red", "blue", "black"), 
	cex = 1, lwd=c(2,1,2,1,1))

dev.off()

